Title

Data

Parameters

suffix = "01_22"
data_to_read = ""

functions

source_from_github(repositoy = "DEG_functions",version = "0.2.1")
source_from_github(repositoy = "HMSC_functions",version = "0.1.0",script_name = "functions.R")

MYB expression

all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)

CNV

new.cluster.ids <- c("cancer", #0
                     "cancer", #1
                     "CAF", #2
                     "cancer", #3
                     "Endothelial", #4
                     "cancer", #5
                     "cancer", #6
                     "CAF", #7
                     "CAF", #8
                     "CAF", #9
                     "cancer", #10
                     "CAF", #11
                     "cancer", #12
                     "cancer", #13
                     "cancer", #14
                     "cancer", #15
                     "cancer", #16
                     "WBC", #17
                     "CAF" #18
                     )

#rename idents:
acc_all_cells = SetIdent(object = acc_all_cells,value = "RNA_snn_res.1") #active snn graph
names(new.cluster.ids) <- levels(acc_all_cells) #add snn graph levels to new.cluster.ids
acc_all_cells@meta.data[["seurat_clusters"]] = acc_all_cells@meta.data[["RNA_snn_res.1"]]
acc_all_cells = SetIdent(object = acc_all_cells,value = "seurat_clusters")
acc_all_cells <- RenameIdents(acc_all_cells, new.cluster.ids) 

# divide "cancer" into patients:
cell_types = acc_all_cells@active.ident %>% as.data.frame()
cell_types[,1]<- as.character(cell_types[,1])
cell_types = cbind(cell_types,acc_all_cells$patient.ident) %>% setnames(old = names(.), 
         new = c('cell_type','patient'))
cell_types[cell_types$cell_type == "cancer",] = cell_types[cell_types$cell_type == "cancer",2]


# hmsc_rows = (startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# acc_rows = !(startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# cell_types[,1][hmsc_rows]  = "HMSC"
# cell_types[,1][acc_rows]  = "ACC"

#add to metadata:
cell_types[,2] = NULL 
cell_types[cell_types$cell_type == "ACC1",] = "HMSC"
acc_all_cells = AddMetaData(object =acc_all_cells ,metadata = cell_types,col.name = "cell.type")

CNV plot

Original score

original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )

myoscore=apply(all_acc_cancer_cells@assays[["RNA"]][original_myo_genes,],2,mean)
lescore=apply(all_acc_cancer_cells@assays[["RNA"]][original_lum_genes,],2,mean)
all_acc_cancer_cells=AddMetaData(all_acc_cancer_cells,lescore-myoscore,"luminal_over_myo")
FeaturePlot(object = all_acc_cancer_cells,features = "luminal_over_myo")

data = FetchData(object = all_acc_cancer_cells,vars = "luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=luminal_over_myo)) + 
  geom_density() 
)

myoscore=apply(acc1_cancer_cells@assays[["RNA"]][original_myo_genes,],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]][original_lum_genes,],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "luminal_over_myo")
print(
    data %>% 
    ggplot(aes( x=luminal_over_myo)) + 
    geom_density() 
  )

0.35 Most correlated score

Myo genes

myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.35)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  20"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "COL16A1"     "RP1-39G22.4" "ACTG2"       "CD200"       "MYLK"        "TP63"        "KCNMB1"      "ADAMTS2"     "LOXL2"      
[10] "TPM2"        "CLIC3"       "SNCG"        "ACTA2"       "TAGLN"       "A2M"         "NGFR"        "CNN1"        "PPP1R14A"   
[19] "MYL9"        "POM121L9P"  

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

intersect(top_myo,original_myo_genes) 

[1] "MYLK"  "TP63"  "ACTA2" "TAGLN"

myo_enrich_res = genes_vec_enrichment(genes = top_myo,background = rownames(acc1_cancer_cells),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)

myo_enrich_res

Lum genes

lum_protein_markers = c("KIT")
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.35)

Warning in cor(expression %>% t(), markers_average) :
  the standard deviation is zero

print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  22"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

 [1] "FBXO44"        "USP48"         "MPC2"          "SLC19A2"       "ATL2"          "B3GNT2"        "AC093162.5"   
 [8] "MAP2"          "KIT"           "GLRB"          "EFNA5"         "PCDHGB7"       "SLC29A1"       "SASH1"        
[15] "ALDH3B2"       "CCND1"         "RP11-254B13.1" "VSIG10"        "NDFIP2"        "SUSD6"         "SPPL2A"       
[22] "DSC2"         

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

intersect(top_lum,original_lum_genes) 

[1] "KIT"

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = rownames(acc1_cancer_cells),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

top correlated score

myoscore=apply(acc1_cancer_cells@assays[["RNA"]][top_myo,],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]][top_lum,],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"top_cor_luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "top_cor_luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "top_cor_luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=top_cor_luminal_over_myo)) + 
  geom_density() 
)

enriched genes score

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

myoscore=apply(acc1_cancer_cells@assays[["RNA"]]@data[myo_enriched_genes,],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]]@data[lum_enriched_genes,],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"enriched_luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "enriched_luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "enriched_luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=enriched_luminal_over_myo)) + 
  geom_density() 
)

0.2 Most correlated score

myo Genes

myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.2)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  473"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "MIR429"        "EXOSC10"       "PLOD1"         "TNFRSF8"       "PDPN"          "FBLIM1"        "RP11-276H7.3" 
 [8] "HSPG2"         "EPHA8"         "IFNLR1"        "GRHL3"         "RP3-426I6.6"   "SDC3"          "RP11-490K7.1" 
[15] "RP11-73M7.1"   "TINAGL1"       "COL16A1"       "RP1-117O3.2"   "RP11-435D7.3"  "RP1-39G22.4"   "SCMH1"        
[22] "RP5-850O15.4"  "RP5-866L20.3"  "NEXN"          "RP4-714D9.2"   "RP5-837M10.3"  "RP5-964H19.1"  "CSF1"         
[29] "RP11-498A13.1" "RP11-96K19.2" 

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

intersect(top_myo,original_myo_genes) 

[1] "MYLK"  "TP63"  "ACTA2" "DKK3"  "TAGLN" "CDH3" 

myo_enrich_res = genes_vec_enrichment(genes = top_myo,background = rownames(acc1_cancer_cells),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)

myo_enrich_res

Lum Genes

lum_protein_markers = c("KIT")
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.2)

Warning in cor(expression %>% t(), markers_average) :
  the standard deviation is zero

print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  1535"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

 [1] "RP11-206L10.11" "HES4"           "RP11-54O7.18"   "CPTP"           "RP3-395M20.8"   "TNFRSF14"       "RPL22"         
 [8] "ICMT"           "ERRFI1"         "RERE"           "H6PD"           "TMEM201"        "NMNAT1"         "UBE4B"         
[15] "APITD1-CORT"    "CENPS"          "RP4-635E18.7"   "FBXO2"          "FBXO44"         "UQCRHL"         "MST1P2"        
[22] "SDHB"           "PADI2"          "ARHGEF10L"      "AKR7A2"         "TMCO4"          "EIF4G3"         "USP48"         
[29] "NIPAL3"         "PIGV"          

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

intersect(top_lum,original_lum_genes) 

[1] "KIT"    "EHF"    "LGALS3"

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = rownames(acc1_cancer_cells),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

correlated genes in original score

myo_intersected = intersect(top_myo,original_myo_genes) 
lum_intersected = intersect(top_lum,original_lum_genes) 
message("genes in myo score:")

genes in myo score:

myo_intersected

[1] "MYLK"  "TP63"  "ACTA2" "DKK3"  "TAGLN" "CDH3" 

message("genes in lum score:")

genes in lum score:

lum_intersected

[1] "KIT"    "EHF"    "LGALS3"

myoscore=apply(acc1_cancer_cells@assays[["RNA"]]@data[c("MYLK","TP63" ,"ACTA2", "DKK3", "TAGLN", "CDH3" ),],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]]@data[c("KIT" ,"EHF",  "LGALS3"),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"enriched_luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "enriched_luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "enriched_luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=enriched_luminal_over_myo)) + 
  geom_density() 
)

enriched genes

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

 [1] "FBLIM1"  "NEXN"    "NMNAT2"  "MYLK"    "CCDC50"  "IGFBP7"  "RAI14"   "ARAP3"   "SPARC"   "CALD1"   "LOXL2"   "COL5A1" 
[13] "ACTA2"   "DKK3"    "MSRB3"   "COL4A1"  "ACTN1"   "TPM1"    "TGFB1I1" "ADORA2B" "MXRA7"   "CNN1"    "TP63"    "ACTA2"  

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

 [1] "PADI2"      "PATJ"       "PEX11B"     "APH1A"      "C1orf43"    "EFNA4"      "NECTIN4"    "SCYL3"      "ELF3"      
[10] "SOX13"      "IRF6"       "MED28"      "EPB41L4A"   "RGL2"       "C6orf132"   "TPD52L1"    "ICA1"       "MACC1"     
[19] "TRPS1"      "FAM83H"     "RASEF"      "ARRDC1"     "COMMD3"     "ANK3"       "GSTO2"      "PDCD4"      "EHF"       
[28] "ALDH3B2"    "SHANK2"     "SORL1"      "FKBP4"      "PTPN6"      "DUSP16"     "RERG"       "ADCY6"      "ERBB3"     
[37] "ERP29"      "SUSD6"      "RPS6KA5"    "SPINT1"     "FEM1B"      "TLE3"       "SCAMP2"     "CLN3"       "ADGRG1"    
[46] "ATP2C2"     "GGT6"       "MYO1D"      "ST6GALNAC2" "CYB5A"      "BLVRB"      "VRK3"       "SYCP2"      "TMPRSS2"   
[55] "LIMK2"      "KIT"       

myoscore=apply(acc1_cancer_cells@assays[["RNA"]]@data[myo_enriched_genes,],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]]@data[lum_enriched_genes,],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"enriched_luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "enriched_luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "enriched_luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=enriched_luminal_over_myo)) + 
  geom_density() 
)

enriched genes and in original score

myo_enriched_genes = myo_enriched_genes[myo_enriched_genes %in% original_myo_genes]
lum_enriched_genes = lum_enriched_genes[lum_enriched_genes %in% original_lum_genes]

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

[1] "MYLK"  "ACTA2" "DKK3"  "TP63"  "ACTA2"

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

[1] "EHF" "KIT"

myoscore=apply(acc1_cancer_cells@assays[["RNA"]]@data[myo_enriched_genes,],2,mean)
lescore=apply(acc1_cancer_cells@assays[["RNA"]]@data[lum_enriched_genes,],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,lescore-myoscore,"enriched_luminal_over_myo")
FeaturePlot(object = acc1_cancer_cells,features = "enriched_luminal_over_myo")

data = FetchData(object = acc1_cancer_cells,vars = "enriched_luminal_over_myo")
print(
  data %>% 
  ggplot(aes( x=enriched_luminal_over_myo)) + 
  geom_density() 
)

HPV

Only HMSC cancer cells:

data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
print(
  data %>% 
  ggplot(aes( x=HPV33.reads)) + 
  geom_density()
)